Mathematica memo

basics

Ctrl + Q で終了
SHIFT + ENTER で実行
Directory[] : = pwd
花文字のI Esc scI Esc
ラグラジアン Esc scL Esc
チルダ Ctrl &, then ~ to enter directly
下付きと上付き両方：下付きを入力後に %

関数のオプションを知りたい。例えば Integrate の optionを知りたいとき
Options[Integrate]

丸めたり四捨五入するとき
Round[] : 四捨五入
Floor[] : くりさげ
Ceiling[] : くりあげ

例
Round[2.7]=3
Round[-2.7]=-3
Floor[2.7]=2
Floor[-2.7]=-3
Ceiling[2.7]=3
Ceiling[-2.7]=-2

Show

LogPlot をいくつかしたあとに Show でグラフを重ねるとする。そのとき


PlotRange ->{All, {10^(-10), 10^(-5) }}

とすると表示がおかしい。次のようにすればよい。


PlotRange ->{All, {Log[10^(-10)], Log[10^(-5)] }}

データ読み込み

"hoge.dat" に次のように書かれていたとする。

#mDM[GeV]  coupling   xsec    omegah2  c1   c2
100        0.1        100     0.120     0     0
200        0.2        101     0.121     0     0
300        0.3        102     0.121     0     0
400        0.4        103     0.122     0     0
500        0.5        104     0.122     0     0

これを読み込むにはImport["hoge.dat","Table"]とすればよい。

読み込みに使うmathematica file と同じディレクトリに dat ファイルを置いている場合、

SetDirectory[NotebookDirectory[]]
data=Import["hoge.dat","Table"];

とすればよい。

違うところにある場合はフルパス書けば間違いない。例えば、

data=Import["/Users/abetomo/work/hoge.dat","Table"];

これで data には

{{#mDM[GeV], coupling, xsec, omegah2, c1, c2},
{100, 0.1, 100, 0.120, 0, 0},
{200, 0.2, 101, 0.121, 0, 0},
{300, 0.3, 102, 0.121, 0, 0},
{400, 0.4, 103, 0.122, 0, 0},
{500, 0.5, 104, 0.122, 0, 0}}

と入る。しかし、１行目は、この後色々やる際に邪魔である。次のようにして消しておく。

data=Delete[data,1]

６つパラメータがあるが、２次元プロットするなら２つだけ欲しい。例えば mDM vs xsec の図を作りたければ、１列目と３列目をとってくれば良い。それは次のようにやるのが速い。(for文でもできるが遅いのでお勧めしない)

data1=Transpose[{Transpose[data][[1]],Transpose[data][[3]]}]

[[i]] は i 行目をとってくるので、Transposeで転置をとって[[i]] すればもとの i列目をとってこれる。最後に Transpose し直す。 Transpose[A] は A Esc t Esc と打つと楽かもしれない。あとは ListPlot[data1] でmDM vs xsec の図ができる。
mDM と coupling の関数 f(mDM, coupling) を mDM の関数としてプロットしたいとする。 f(x,y) はどこかで定義しておく。そのときはMapThreadを使えばよい。

In[1]:= MapThread[f,{{a1,a2,a3},{b1,b2,b3}}]
Out[1]= {f[a1,b1], f[a2,b2], f[a3,b3]}

これを利用し、次のようにすれば、mDM vs f(mDM, coup) の絵が描ける。

f[x_,y_]:=x+y
mDMtab = Transpose[data][[1]];
coupltab = Transpose[data][[2]];
ftab=MapThread[f,{mDMtab,coupltab}];

dam=Transpose[{mDMtab, ftab}];
ListPlot[dam]

interpolation

{
 {{x0,y0},z0},
 {{x1,y1},z1},
 {{x2,y2},z2},
 ...
}

f = Interpolation[%]

f[x,y]

interpolation するファイルの作り方

What we want is

{
 {{x0,y0},z0},
 {{x1,y1},z1},
 {{x2,y2},z2},
 ...
}

The 1st method

dam1={x0, x1, x2,...}
dam2={y0, y1, y2,...}
dam3={z0, z1, z2,...}

dam4 = Transpose[{dam1, dam2}]
Table[{dam4[[i]], dam3[[i]]}, {i, 1, Dimensions[dam3][[1]]}]

The 2st method

dam1={x0, x1, x2,...}
dam2={y0, y1, y2,...}
dam3={z0, z1, z2,...}

Transpose[{Transpose[{dam1, dam2}],dam3}]

グラフを並べる

GraphicsGrid[
 Table[
    DensityPlot[Sin[x + a] Cos[y + b], {x, 0, Pi}, {y, 0, Pi}]
   ,{a, 0, 2}, {b, 0, 2}
    ]
]

GraphicsGrid[{{plot1,plot2}}]

Plot の Plotlabel を右上に表示させる

デフォルトだと図の真上に表示される。横軸の名前だと思われたことがあるので、右はしに寄せたい。 Pane をつかう。

   RegionPlot[x^2 + y^2 < 1
    , {x, -1, 1}, {y, -1, 1}
    , PlotLabel -> Pane["This is the title", Alignment -> Right, ImageSize -> 300]
    , ImageSize -> 300 
   ]

ここで、ImageSize を Pane と Plot の両方で指定しないといけない。

データファイルでゴミを落とす

例えば、 dam = {{1,2,3},{1,2,5},{1,3,a}} とあったら、３行目は a が数字が入ってないので、これを落としたいときがある。このときは、

max = Dimensions[dam][[1]];
damdel = {}
For[i = 1, i < max + 1, i++,
 If[!NumberQ[dam[[i, 3]]],
  Print[i];
  damdel = Append[damdel, {i}]
  ]
 ]

これで、damdelには数字が入っていない行番号が入る。今の例だと、{{3}}　となる。たくさんある場合には、 {{3},{6},{10},...} などのようになる。あとは

damnew=Delete[dam, damdel]
Clear[max, damdel,i]

とすればdamnewには数字が入っていないところを消したものがはいる。

データを出力ときに簡単にする

数字をデータファイルに出力するとき、有効数字３桁くらいで欲しいのに、10桁くらい出力されて非常に見づらい。この場合は、fortran や C で書くように、桁数を制御したい。以下のようにすれば何とかなる。 reference

f77Eform[x_?NumericQ, fw_Integer, ndig_Integer] := 
  Module[{sig, s, p, ps}, {s, p} = MantissaExponent[x];
   {sig, ps} = {ToString[Round[10^ndig Abs[s]]], ToString[Abs[p]]};
   StringJoin @@ 
    Join[Table[" ", {fw - ndig - 7}], If[x < 0, "-", " "], 
     If[StringLength[sig] < ndig + 1, {{"0."}, {sig}, 
       Table["0", {ndig - StringLength[sig]}]}, {{"1."}, 
       Table["0", {StringLength[sig] - 1}]}], {"e"}, 
     If[p < 0, {"-"}, {"+"}], 
     Table["0", {2 - StringLength[ps]}], {ps}]];

f77[x_?NumericQ] := f77Eform[x, 10, 3]; (<---- １０文字分確保して、有効数字は３桁)

Export["st_pre.dat", f77 /@ {sparapre,tparapre,contourpre}, "Table"];

半透明なグラフのpdf出力にメッシュを出さない (1)

png で出力すれば解決するが、やはりベクター形式のきれいな図が欲しい。そのためには面倒であるが以下のステップを踏む必要がある。まず、図を作るときに、

   ContourPlot[f[x, y]
    , {x, 10, 350}, {y, 10, 3500}
    , ContourShading -> {{Orange, Opacity[0.5]}, None}
    , MeshStyle -> None
   ]

とやると、Mathematica 上ではメッシュは出ない。半透明にするときは、他の図と重ねる場合なので show コマンドを使うが、

Show[..., Method -> {"TransparentPolygonMesh" -> True}]

とする必要がある。するとメッシュ無しの図が画面に出力される。

これを pdf にする場合、右クリックして画像を保存を選んではいけない。これをやるとメッシュが現れる。メッシュを出力しないためには、右クリックして Print Graphic... を選び、そこから pdf に出力する。

この方法だとメッシュは出ないが、大きな余白がついてくる。余白を消すには、ターミナルで "pdfcrop hoge.pdf" とする。これで余白の無い "hoge-crop.pdf" が手に入る。 Mac で Mathematica 9 の場合なので、環境によってはこんなこと必要ないかもしれない。

半透明なグラフのpdf出力にメッシュを出さない (2)

上記の方法は面倒くさい。いちいち右クリックしないといけないし、さらにコマンドを打たないといけない。これだと大量に図を処理するときに大変である。

一発でやれる方法を見つけた（参考）。
図を生成する前に以下のコマンドを実行しておくだけである。

Map[SetOptions[#, 
    Prolog -> {{EdgeForm[], Texture[{{{0, 0, 0, 0}}}], 
       Polygon[#, VertexTextureCoordinates -> #] &[{{0, 0}, {1, 
          0}, {1, 1}}]}}] &, {Graphics3D, ContourPlot3D, 
   ListContourPlot3D, ListPlot3D, Plot3D, ListSurfacePlot3D, 
   ListVectorPlot3D, ParametricPlot3D, RegionPlot3D, RevolutionPlot3D,
    SphericalPlot3D, VectorPlot3D, RegionPlot, ContourPlot}];

あとは普通に図を生成して Export["hoge.pdf", %] とすれば良い。 Mesh->None とかは必要ないようだ。

RegionPlot が挙動不審

RegionPlot で使う関数の中で NSolve を使っていたら動かなかった。Mathematica9 では動いたが Mathematica12 では動かない。その場合はRegionPlotのオプションで以下を試す。

  RegionPlot[ ....
  , "NumericalFunction" -> False
  ]

ContourPlot とかでもおかしいときには試すと動く場合がある。

Log plot in Contour Plot

reference
For the case the linked page is erased, we copy the information there below.

==================================================
You need to do the scaling by hand, using FrameTicks. Here is a simple 
example:

f[x_,y_]:=x^2+Log[y]^2

Let's plot this over the range -3< x <3 and Exp[-3]< y < Exp[3] so that we 
get a nice circle.

xmin=-3;
xmax=3;
ymin=Exp[-3];
ymax=Exp[3];

Now, we want to use ContourPlot with Log scaling. So

ContourPlot[f[x, 10^z],{x,xmin,xmax},{z,Log[10,ymin],Log[10,ymax]}]

will produce the desired plot, but with ticks labeled using the log 
scale. We fix the ticks by using FrameTicks. Here is a function that 
produces a table of {linear,log} values:

tickfunction[min_,max_]:=
Flatten[
   N@Table[
     {e+Log[10,d],d 10^e},
     {e,Floor[Log[10,min]],Ceiling[Log[10,max]]},
     {d,1,5,4}
   ],
   1
]

So, the following produces the desired plot:

ContourPlot[f[x,10^z],{x,xmin,xmax},{z,Log[10,ymin],Log[10,ymax]},
    FrameTicks->{Automatic,tickfunction[ymin,ymax]}]

Carl Woll
Wolfram Research
==================================================

The above function gives like {0, 0.5, 1, 5, 10, 50, ...}. To make it like {0.1, 1, 10, 1000, ...}, use the following.

tickfunction[min_, max_] :=
 Flatten[
  Table[
   {N[e + Log[10, d]], 
    If[d == 1, Round[d 10^e], ""], {If[d == 10, 0.018, 0.009], 0}},
   {e, Floor[Log[10, min]], Ceiling[Log[10, max]]},
   {d, 1, 10, 1}
   ],
  1
  ]

tickfunction2[min_, max_] :=
 Flatten[
  N@Table[
    {e + Log[10, d], "", {If[d == 1, 0.018, 0.009], 0}},
    {e, Floor[Log[10, min]], Ceiling[Log[10, max]]},
    {d, 1, 10, 1}
    ],
  1
  ]

ContourPlot[
 f[10^z, y]
 , {z, Log[10, ma0min], Log[10, ma0max]}
 , {y, 0, 3} 
 , FrameTicks -> {
     {Automatic, Automatic}
    ,{tickfunction[ma0min, ma0max], tickfunction2[ma0min, ma0max]}
   }
 ]

Round[d 10^e] is used for "10. --> 10". If you want to show numbers smaller than 1, modify it.
0.018 and 0.009 is the length of the major and minor ticks.
Note that FrameTicks is

FrameTicks -> {{left, right},{top, bottom}}